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A FRACTIONAL LAPLACE EQUATION: REGULARITY OF 
SOLUTIONS AND FINITE ELEMENT APPROXIMATIONS 

GABRIEL ACOSTA AND JUAN PABLO BORTHAGARAY t 


Abstract. This paper deals with the integral version of the Dirichlet homogeneous fractional 
Laplace equation. For this problem weighted and fractional Sobolev a priori estimates are provided in 
terms of the Holder regularity of the data. By relying on these results, optimal order of convergence 
for the standard linear finite element method is proved for quasi-uniform as well as graded meshes. 
Some numerical examples are given showing results in agreement with the theoretical predictions. 
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1. Introduction. In the last years, the study of nonlocal operators has been 
an active area of research in different branches of mathematics, and these operators 
have been employed to model problems in which different length scales are involved. 
Anomalous diffusion phenomena are ubiquitous in nature [29, 35]; among the several 
applications of these nonlocal models let us mention image processing [8, 22, 24, 32], 
finance [11, 13], electromagnetic fluids [34], peridynamics [42] and porous media flow 
[2, 14]. 

In this work we will be interested in the fractional Laplace operator of order s, 
which we will denote by (—A) s and simply call the fractional Laplacian. In the theory 
of stochastic processes, this operator appears as the infinitesimal generator of a stable 
Levy process [3, 45]. 

If the domain under consideration is the whole space K ra , then there is a natural 
way to define it as a pseudodifferential operator of symbol |£| 2s . Indeed, for a function 
u in the Schwartz class S, let 

(1.1) (-AYu = T- 1 (|£| 2s JFu) , 

where J- denotes the Fourier transform. It is expectable for this operator to approxi¬ 
mate the usual Laplacian for s —> 1 and the identity as s —> 0. 

The fractional Laplacian can equivalently be defined by means of the identity [16] 


where the normalization constant 


C(n, s ) 


2 2 s sr(s + f) 
W 2 r(l — s) 


is taken in order to be consistent with definition (1.1). 

One of the main difficulties arising in the study of this operator is its nonlocality; 
in order to localize it, Caffarelli and Silvestre [9] showed that it can be realized as a 
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Dirichlet-to-Neumann operator by means of an extension problem in the half-space 



However, there is more than one way to define the fractional Laplacian on an 
open bounded set Q: 

1. One first possibility is to consider fractional powers of the Dirichlet Laplace 

operator in the sense of spectral theory. Indeed, let {i/jk, ^k}k&n C x 

K + be the set of normalized eigenfunctions and eigenvalues for the Laplace op¬ 
erator in fl with homogeneous Dirichlet boundary conditions, so that {il>k}keN 
is an orthonormal basis of L 2 (H) and 

f = Afc'i/’fc in 

\ ipk = 0 on dtt. 

Then, the spectral fractional Laplacian (— A) a s is defined for u £ C£°(Q) by 

OO 

(-A) s s u = ^2{u,'ip k )\ s k u k , 

k =l 

and can be subsequently extended by density to the Hilbert space H S (H). In 
[44], the Caffarelli-Silvestre result was proved for this operator, thus achieving 
a local problem posed on a semi-infinite cylinder Q x (0, oo). This localization 
was exploited by Nochetto, Otarola and Salgado in [37], where the authors 
study the numerical approximation of the spectral fractional Laplacian by 
considering graded meshes in the extended variable. 

2. A second feasible definition is attained by considering the integral formula¬ 
tion (1.2), and restricting it to functions supported in H. This gives rise to 
the integral fractional Laplacian (—A)|w. This operator is different to the 
spectral fractional Laplacian; for example, their difference is positive definite 
and positivity preserving [36]. See also [41], where the spectra of these op¬ 
erators are compared. The main difficulties to overcome when dealing with 
numerical analysis of this integral fractional Laplacian are associated to its 
nonlocality and to the singularity at x = y of the kernel it involves. 

3. Finally, it is also possible to consider a regional fractional Laplacian (— A) S R , 
in which the integration in (1.2) is restricted to H. This operator is known to 
be the infinitesimal generator of the so-called censored stable Levy processes 

[4]- 

Throughout this work, we will restrict to the second definition, and for sake 
of simplicity we will skip the I subindex when denoting it. The analysis we will 
perform is valid for the regional fractional Laplacian as well. Moreover, it can be 
straightforwardly extended to operators of the form 

Cjcu(x) = P.V. [ ( u(x ) - u(y))K(x - y) dy , 

Jr™ 

where the kernel K : R" \ {0} —> (0, oo) satisfies 

jK £ L 1 (R n ), with j(x) = min{ 1, |x| 2 }, 

39 > 0, s £ (0,1) such that K{x) > 0\x\~^ n+2s \ ieM" \ {0}, 

K(x) = K(—x), i£M" \ {0}. 
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Numerical approximation for the fractional Laplacian on bounded domains has 
been addressed in the last years. D’Elia and Gunzburger [15] exploited the nonlocal 
vector calculus introduced in [17] in order to perform a study of convergence of certain 
approximations of the fractional Laplacian as the nonlocal interactions become infi¬ 
nite. Huang and Oberman [27] proposed a method which combines finite differences 
with numerical quadrature, obtained a discrete convolution operator and studied nu¬ 
merically the convergence and order of their method in the L°°(Cl) norm. However, 
these methods were implemented only in 1 dimension, and regularity of solutions in 
those works is assumed as part of the hypotheses. 

This work is an attempt to deal with basic analytical aspects required to convey 
a complete finite element analysis of the following fractional Laplace problem 


(1.3) 


f (—A ) s u = f in Cl, 
{ zt = 0 in Cl c , 


where 0 < s < 1 and (—A ) s u denotes the operator defined in (1.2). Our aim is 
to provide regularity of solutions as well as a priori error estimates for the discrete 
approximations together with a feasible finite element implementation. 

From now on we assume that Cl C R” is a bounded Lipschitz domain and / a 
bounded function defined in Cl. 

The fractional Laplace equation in the form (1.3) shares some key analytical 
features of the classical Laplacian making it amenable - in principle - to a direct finite 
element treatment. Nevertheless, Sobolev regularity results for this problem are recent 
and expressed in terms of Hormander spaces [26]. Moreover, in that paper the domain 
is required to have a C°° boundary, which makes its results not entirely satisfactory 
for a FE analysis. Also, some numerical difficulties -such as handling the singularity 
of the kernel- seem to be the main disadvantages that have discouraged a direct FE 
approach. Concerning the latter, we found that applying rather standard techniques 
(actually borrowed from the theory of the Boundary Element Method [39]) together 
with an appropriate treatment of the integrals involving the unbounded domain Cl c , 
accurate FEM solutions (at least in 2D) could be delivered. 

In this paper we address Sobolev regularity results for (1.3) in less regular do¬ 
mains. We provide weighted fractional Sobolev a priori estimates in terms of the 
Holder regularity of the data for Lipschitz domains satisfying the exterior ball condi¬ 
tion. As it is shown in Remark 3.13, our predicted regularity is sharp. The proof relies 
on recent Holder regularity results of Ros-Oton and Serra [38] (see Section 3). Even 
though Sobolev-Sobolev (instead Sobolev-Holder) results are preferable, this theorem 
is a key tool within the FE analysis developed along this work. 

The nonlocal nature of problem (1.3) is reflected in the fact that fractional Sobolev 
norms are not additive respect to the domains. Therefore, after finding a suitable in¬ 
terpolation operator and getting adequate local interpolation estimates, some cautions 
must be taken in order to bound the global interpolation error [19, 20]. In order to 
deal with graded meshes we extend well known error estimates for the Scott-Zhang 
interpolation operator to our weighted fractional Sobolev spaces. These estimates are 
derived by introducing improved Poincare inequalities in the fractional setting in a 
form that we were unable to find in the literature. 

The organization of this paper is as follows. Section 2 is devoted to basic defi¬ 
nitions and properties of the operator and spaces involved. In Section 3, the matter 
of regularity of solutions is addressed; it is proved that under some Holder regularity 
hypotheses on the source /, the solutions of the problem under consideration gain half 
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a derivative in the Sobolev sense. Moreover, if the source is more regular, sharper 
estimates in weighted fractional spaces are shown to hold. Then, in Section 4 interpo¬ 
lation estimates are developed both in standard and weighted spaces, thus obtaining 
a priori error estimates for Finite Element (FE) solutions of equation (1.3). Finally, 
in Section 5, some remarks on the implementation are discussed and several numer¬ 
ical tests are presented showing results in complete agreement with our theoretical 
predictions. More details about the implementation in dimension two can be found 
in the forthcoming paper [1], and related results for the fractional eigenvalue problem 
in [5], 

Acknowledgments. The authors would like to thank L. Del Pezzo and S. 
Martinez for stimulating discussions on the topic and valuable help by thoroughly 
reading draft versions of this paper. 

2. Analytic setting. In this section we set the notation and review some prop¬ 
erties of the spaces involved in the rest of the paper. Let us begin by recalling some 
function spaces. Throughout this work, s is a parameter such that 0 < s < 1 and 
C K™ (n > 1) is a bounded Lipschitz domain. 

2.1. Function spaces. The Sobolev space H s (fl) is defined by 
H s (n) = {»£ L 2 (D): |u| ff , ( n) < oc} , 


where 





(u(x) - v(y)) 2 
\x - y\ n+2s 



2 


denotes the Aronszajn-Slobodeckij seminorm. The space H S (Q) is a Hilbert space, 
equipped with the norm 


IMI-R s (n) — IMI,L 2 (n) + Mu^n)- 

Equivalently, the space H S (Q) may be regarded as the restriction to f l of functions in 
H s (M. n ). Zero trace spaces H^(n) can be defined as the closure of w.r.t. the 

H s norms. Equivalently, if the boundary of the domain is smooth, they can be defined 
through real interpolation of spaces by the A-methocl (for example, [31, Chapter 1]). 
It is well-known that for 0 < s < 1/2, the identity holds. 

Sobolev spaces of order grater than 1 are defined in the following way: given 
k £ N, then 


H k+S (ty = {v e H k (n): \D a v\ € H s (n)Va s.t. |a| = k} , 
furnished with the norm 

IHIffM-s(fi) = IMIfpqn) + ^2 \D a v\ H s(ny 

\at\=k 

In the sequel we assume that k = 0 or k = 1, which are the cases of interest along our 
presentation. 

Let us recall that weighted Sobolev spaces are a customary tool for dealing with 
singular solutions. In the present context we find useful some fractional and weighted 
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spaces. Our weights are powers of the distance to the boundary of Cl. Accordingly, 
we introduce the notation 6 : R" —> R> 0 , S(x) = d(x,dCl) and consider the norm 


HI 


Hi 


s (n) 



\D^v(x) — D^v(y) | 2 

X — y\ n + 2s 


S(x, y) 2a dx dy 


where a > 0, s G (0,1) and 


S(x, y) = min{(S(a;), <%)}. 


In this way we write 


Hi(Q) = {v € H\Q): |M| ff « (n) < oo} . 

Global versions H l a n (R") are easily obtained integrating in the whole space R" and 
taking <5 as before. In the sequel we drop the reference to Cl in the global case and 
write H e a (R n ). 

Remark 2.1. Although we are interested in the case a > 0, let us recall that in 
the definition of standard weighted Sobolev spaces H*( O), with k being a nonnegative 
integer, arbitary powers of 5(x) can be considered [30]. On the other hand, for general 
weights some restrictions must be taken into account in order to get a right definition 
of the spaces. A classical family of weights is that of the Muckenhoupt A 2 class [f3[. 
In the global version H^(R n ) we restrict the range of a to 0 < a < 1/2 in order 
to have S a G A 2 . Moreover, this restriction arises naturally in our results involving 
graded meshes (see Remark 5.2). 

Weak solutions of equation (1.3) may be defined by multiplying by a test function 
and integrating by parts the Laplacian term. Namely, consider the space 

Y = {ueH s ( IT) : u = 0 in fi c }, 


equipped with the H s ( 1") norm. For V it is known that C'o°(f2) is a dense set [21]. 
Then, the weak formulation of (1.3) is: find u G V such that 


( 2 . 1 ) 


C(n,s) ff (u(x) - u(y))(v(x) - v(y)) 


- y I 


n+2s 


dxdy= / f(x)v(x)dx 
Jn 


for all v G V, where Q = (fl x R") U (R n x f l). The integral in the right hand side 
of the previous equation makes sense if / belongs to the dual space of Y, which we 
denote by V*. Observe that (iP(R"))* C V*. 

At this point we recall for further reference the following result [7]. 

PROPOSITION 2.2 (Poincare inequality I). Let S be an star-shaped domain w.r.t. 
a ball B. For any u G H S (S) with 0 < s < 1 we call v = f s v. Then we have 


Ik-v|H(S) < cd|k|H-» ( S), 


with c bounded in terms of (g), where ds = diam(S) and ds = diam(B). 
Proof. We write 





(v( x ) 


v{y))dy 


2 

dx < 



v(y)) 2 dydx, 
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therefore 


L 


(v — v) 2 dx < 


d n s +2s 

|S| 


{y(x)~ v{y)) J J 
\ x -y\n + 2 s dydx- 


Taking into account that dg ~ \B\ < |S|, the Proposition follows. □ 

Remark 2.3. Following [7] we call the chunkiness parameter of S. Another 
well-known result is the following. We include its proof here for completeness. 

Proposition 2.4 (Poincare inequality II). Given LI as above, there exists a con¬ 
stant c = c(Ll, n , s) such that 


IMU 2 (n) < cMtfspun) VveV 


Proof. By Lemma 6.1 of [16], there exists some constant c(n, s) > 0 such that for 
all x £ n, 


c(n,s)|0| » < f - -- 

Jn° \x-y 


n-\-2s 


dy. 


On the other hand, since v = 0 on Ll c we know that v(x) 2 = (v(x) — v(y)) 2 for all 
x e LI, y £ f2 c . So, we obtain 


z{n, s)|f2| " f v{x) 2 dx< ff 
J n JJ n 


{v{x) - v(y)f 


QxO c \ x 


_ y\n+2s 


dx dy, 


and the Poincare inequality follows straightforwardly. □ 

An immediate consequence of Proposition 2.4 is that the bilinear form a : Yx V —> 
R, 


( 2 . 2 ) 


a(u, v) = 


C(n,s ) ff (u(x) - u(y))(v(x) - v(y)) 


l 


■ - y I 


n+2s 


dx dy 


is coercive. Its continuity is an obvious consequence of the Cauchy-Schwarz inequality. 
Therefore, it can be easily seen, by applying the Lax-Milgram theorem, that if / € Y* 
then there exists a unique u e V solution of problem (2.1). 

As the H s (M n ) seminorm is equivalent to the H s ( R") norm on V, let us define 


II II / \ — (^b ^0 I I 

IMIv := a(v,v)* = y—-— 

The approach we shall follow to obtain error estimates is simply to consider an 
adequate interpolator in a FE space V^, and estimate the interpolation error. In order 
to achieve this, it is convenient to understand the relation between the norm in V, 
and the (semi)norm in H s (Ll). 

PROPOSITION 2.5 (Hardy-type inequalities, see [25, 18]). Let S2 be as above. If 
0 < s < then there exists c = c(f2, n, s) > 0 such that 

(2.3) j dx < c\\u\\ 2 HB{n) Vue H'(Q). 

If ij < s < 1, then there exists c = c(S2, n, s) > 0 such that 

dx < c\u\ 2 H s(^) Vue Co(fl). 


f \u(x )\ 2 

n S(x) 2s 
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As a consequence of the previous Proposition we get the following 
Corollary 2.6. If 0 < s < ^, then there exists a constant c = c(f2, n, s) > 0 
such that 

IMIv<c|M|*. ( n) VreV. 

On the other hand, if i < s < 1 there exists a constant c = c(Cl, n, s) > 0 such that 

IMIv < cMffs(n) v«e V 


Proof. Let us prove the first inequality; the second one follows in the same fashion 
recalling that we can assume v £ Cfi°(f2) [21] and concluding by density arguments. 
Let v € V, then, since v = 0 in Cl c , 


Mlv = 


C(n, s) 




\v{x)-v(y) | 2 , t ff 

\ x _ sr+1 , dxdy + C(n,s)jj^ 


nxn c 


M^)! 2 

\x - y | n+2s 


dxdy 


< c(n,s) \v\ 2 H * {n) + [ \v{x)\ 2 [ 

in jb{x, 


5{x)Y 


\x-y\ 


n+2s 


dy dx 


= c(n,s) \v\ 2 Ha(n) 


I 

in 


Hx )\ 2 

S(x) 2s 


dx 


Applying the Hardy inequality (2.3), the first estimate follows. □ 

3. Sobolev regularity. The purpose of this section is to provide regularity 
estimates for solutions of (1.3) in terms of fractional Sobolev norms. We start by 
reviewing some key results given in [38]. 

Theorem 3.1 (See Prop. 1.1 in [38]). If Cl is a bounded, Lipschitz domain 
satisfyinq the exterior ball condition and f £ L°°(Cl), then any solution u of (1.3) 
belongs to C s (R n ) and 


(3.1) 


IMIc«(R») < C'(^) s )ll/IU~(n)- 


Moreover, if / is Holder continuous, estimates for higher order Holder norms of u 
can be obtained. For 0 < /3, we denote by | • \cfi(n) the C^(Cl) seminorm. For 6 > — /?, 
let us write f3 = k + /?' with k integer and f3' £ (0,1]. Following [38] we define the 
seminorm 


l«= s »p 


x,y^Q 


and the associated norm 



in the following way: for 6 > 0, 


f=0 


IMI/j = H [snpS(xY +9 \D e w(x)\ ) + \w\^ 


W 


while for — /? < 6 < 0, 


Ml 


w 

p 


k 

IMIc- e (fi) + 

e=i 


sup S( x) £+e 

k x£fl 
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It holds, 

Theorem 3.2 (See Prop. 1.4 in [38]). Let Li be a bounded domain and (3 > 0 be 
such that neither /? nor [3 + 2 s is an integer. Let f £ C^(fl) be such that ||/||^ < oo, 
and u £ C S (R”) be a solution of (1.3). Then, u £ C^ +2s (Cl) and 

Nfct<c(a s ,/3)(NI^(^) + II/II? ) )- 

In the next remarks we explore some consequences of the previous theorems written 
in a way useful in the sequel. 

Remark 3.3 (Case 0 < s < |). Taking [3 £ (0,1 — 2s) in Theorem 3.2, we get 
that there exists a constant C(Ll,s,/3 ) such that 

(3.2) ™p o S(x,v)>+- '"W < C (ll/IU-m + ||/||Jf>) - 

Moreover, since (3 < 1, for f £ it is simple to prove that 

||/||J ) <C'(0,s)||/|| c , (n) . 


Remark 3.4 (Case i < s < 1). Considering [3 £ (0,2 —2s), Theorem 3.2 implies 


that 


and 


x,y£fl 


< c(n.,MU 

sup S(x) 1 ~ s \Du(x)\ < C (n,s,/3, Il/H^) . 


In the remainder of this section we show how to use these results to bound Sobolev 
norms of u. In order to do that it is useful to divide OxO into a set in which the 
distance between x and y is bounded below by 8{x,y) and a set in which \x — y | 
is smaller than that. Roughly, for the first set, Holder regularity of the solution is 
enough to control the integrand involved in fractional seminorms of u, as this region 
is away from the diagonal. As for the second one, since the weight involving \x — y\ 
is singular at y = x, some extra term is required in order to control its growth; this 
is obtained by means of Theorem 3.2. 

It is convenient to observe that, given a function v : H x Ll —> K. such that 
v{x,y) = v(y,x ) for all x,y £ Q, the integral of v over SI x O equals 2 times its 
integral over the set 


(3.3) 


A = {( x,y) eSlx!!: d{x,y) = <5 (a;)}. 


We make use of the decomposition mentioned in the previous paragraph by defin¬ 
ing 

(3.4) B = {(x, y) £ A: \x - y\ > (5(x)}. 

Remark 3.5. At this point, let us recall an useful identity regarding integrability 
of powers of the distance to the boundary function. The following holds whenever 
a < 1: 

(3.5) [ 5{x)~ a dx = 0 

.In 

See, for example, the proof of Lemma 2. If in [10]. 
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3.1. Regularity in standard fractional spaces (0 < s < 1/2). We are now 

ready to prove: 

Proposition 3.6. Let 0 < s < | and f £ C^ _s (n). Then, for every e > 0, the 
solution u of (2.1) belongs to H s+ ^~ E (Ll), with 


'H 


,+ i -e (n) — 


C(S2, s, n) 


l c , 7 _s (fi)' 


Proof. Take 9 € (s, 1) and consider the splitting of A mentioned before. Then, 
applying estimate (3.1), 


Kx) - u(y) | ; 
B \x — y\ n + 2e 


dx dy < 


< c(n 


s Jlll llL~(n) 
2 


17 JB(x,S(x)) c 


-y\ 


< 


SJ C 

d-s Jn 


-n—26-\-2s 


dy dx 


A necessary and sufficient condition for the finiteness of the right hand side in 
the previous inequality is that 9 < s + i. 

On the other hand, assume / £ C^(Cl) for some (3 > 0. In a similar fashion the 
application of inequality (3.2) yields 


A\B 




<C <5( x )- 2(/3+s) 


(/ 

\J B(x,S(x )) 


: - y\- n ~ 2e+2l3 + 4s dy dx. 


Now, the integral over B(x,5(x)) is finite if and only if /3 + 2 s > 9. So, in this case 
we obtain 


(3.6) 


I u{x)-u(y)\‘ 
a\b \x - y\ n+2S 


■dxdy <C [ 5(x) 2 ( s e ^dx, 
Jn 


where in the end the constant is of the form 

C(Q,s,n,P) 2 


C = 


i + 2 s-9 


II /II0^(0) ■ 


Once again, the integral in the right hand side of (3.6) is finite if and only if 9 < s + 
If (3 = \ — s, choosing 9 = s + \ — e, we find 



\u(x)-u{y)\ 2 
\x — y \ n + 26 


dxdy < 


C(n,S,n) [ 5(x)~ 1+2e dx. 
£ Jn 


Since the integral in the right hand side is 0(£ _1 ) (recall identity (3.5)), the proof is 
concluded. 0 

Remark 3.7. If f is more regular than C^~ s (Ll), then no further gain of regu¬ 
larity from estimate (3.2) is possible by means of the technique of the previous proof. 
This is indeed sharp, see Remark 3.13. The matter is that the parameter f3 disappears 
in inequality (3.6). 
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3.2. Regularity in fractional weighted spaces (1/2 < s < 1). Next we 
show that an analogous of Proposition 3.6 is possible for 1/2 < s < 1 and hence 
1/2 derivative can be gained in the a priori estimate. Moreover, along the proof of 
this result it becomes clear that the singular behavior of the solution can be localized 
near the boundary. Therefore, introducing appropriate weights we find alternative 
regularity results that can be used to build a priori adapted meshes. 

Before proceeding let us notice that the expected gain of 1/2 derivative would 
imply that the solution belongs at least to 7J 1 (fl). This fact can be proved studying 
the behavior of the fractional seminorms | • |//i-£(n) as e —> 0 . 

Let us recall this useful result proved in [ 6 ]: 

Proposition 3.8. Assume v £ L p ( 12), 1 < p < oo. Then, 

(3-7) lim e|«|^ 1 _., p(n) =C , (n ) p)|«|^ liI1(n) . 


In first place, we want to prove that for the solution u of (2.1), the left hand side 
of (3.7) remains bounded as e —> 0, so that it belongs to if 1 (12). For that purpose, 
we require the following local Holder regularity estimate, (see [38], Lemma 2.9): 
Lemma 3.9. If f £ L°°(12) and 7 £ (0, 2s), then u verifies 


(3-8) M chb**)) ^ CRS ~^\fh-(n) V* € 12, 

where R = and the constant C depends only on 12, s and 7 , and blows up only 
when 7 —> 2s. 

The mentioned H 1 regularity follows from the previous results, and its proof can 
be obtained like the one of Proposition 3.6. 

Lemma 3.10. If 1/2 < s < 1 and f £ L°°(fl), then a solution u of (2.1) belongs 
to H 1 ^) and therefore to iL 1 (R"). Moreover, it satisfies 




C(12, s, r)H/||l°°(o) 
(l-s)(2s-l) 


where the constant C(£l,s,n) is uniformly bounded for all s £ (1/2,1). 

Proof Take s £ (0,1 — s) and in the same fashion as before consider the sets A 
and B, with the slight difference of a instead of a 8{x) in the definition of the 
latter. Taking 7 = 1 — C(e) for some 0 < C(e) < e to be chosen, applying estimate 
(3.8) and proceeding as in the proof of Proposition 3.6, it follows 


\u(x) ~ u(y)\ 2 (7(12,5,71)11/11 

cLy ax ^ ■ 


■(«) 


A\B \x-y\ n + 2 A~e) 


e — C(e) 


[ <5(*) 2(s 

Jn 


' 1+e U: 


X. 


Observe that the constant C in the previous inequality remains bounded for s £ 
( 1 / 2 , 1 ), and that the integral is O (( 2 s — 1 + 2 e) _1 ). 

On the other hand, taking into account the global Holder regularity of u it is 
immediate to obtain 



\u{x)-u{y)\ 2 
\x - y\ n + 2 A~ e ) 


dy dx < 


<7(12,s,n)||/||| 0 


'(n) 


1 — s + e 


[ S(x) 2 ^ s 
Jn 


- 1+£ U: 


X. 


Combining the previous estimates, we get 

■ 1 2 (7(12, s, r)||/|||=o(q) 

H1 ° in> ~ (e — C(e))(l — s + e)(2s — 1 + e) 
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where the constant C(Q,s,n) remains bounded for s £ (1/2,1). Taking C(e) such 
that e — C(e) = 0(e), the desired conclusion follows thanks to Proposition 3.8. □ 

Next, we require some regularity on Du. Let /3 £ (0, 2 — 2s) and assume that 
/ £ C^(fl). Consider the subsets A and B of U x H as before (see eqs. (3.3) and 
(3.4)) and introduce the weighted integral 


I := 



\Du(x) — Du(y)\ 2 
\x - 


8(x,y) 2a dx dy. 


Using the first inequality of Remark 3.4 we explore how to take the involved param¬ 
eters t and a in order to keep I bounded. 


I<C 


f ( f \X- y\W+*s-l)-n-2(l-X) dy \ §{x) 2(o 

J f2 y«/ B(x,S(x)) J 


-P~ s Ux < 


< ___ / S(x) 2{a+s - e) dx < _-_ 

-p + e-2sj n [ ) - (fd + £-2s)(l + 2(a-s-£)y 

where, in order to ensure the convergence of the integrals involved, we must require 

(3.9) £ — (3<2s and l < a + s + 1/2. 


On the other hand, for 

II := 

again due to Remark 3.4, 


| Du(x) — Du(y)\ 2 
B |a; - j,|"+2(<-1) 


S(x,y) a dxdy 


II <C / \x^ y \- n - 2 ^- 1) dy)5{x) 2{a+s - l) dx< 

Jn \JB(x,S(x)) c I 


<C f Sixf^+^dx <-——- 

“ Ja -1 + 2 ( a -s-£)’ 

where the condition for the finiteness of II is guaranteed if we restrict our attention 
to (3.9). Under these conditions, we have proved: 

(3 ' 10) \ Du \Hi(n) < ( /3 + £_ 2s )(i + 2 (a-s-£))' 

Within the range provided in (3.9) we can highlight some cases of interest. In 
the same spirit of Proposition 3.6, we have, considering a = 0 and £ = s + 1/2 — e in 
(3.10): 

Proposition 3.11. 1/1/2 < s < 1 and f £ CP(Q) for some /? > 0 , then the 
solution u of (2.1) belongs to H s+ 2 ~ e (H) for all e > 0, with 

in I ^ C(U, s, n, 0) 

I ~ ^/E(2s-l) 
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If we restrict the weight to the Muckenhoupt A 2 class (see Remark 2.1), which 
can be relevant for extending this considerations to the global case treated later, we 
need to choose a < 1/2. This restriction is also of importance in the optimality of 
the graded meshes proposed later. Accordingly, assume a = 1/2 — e: for e > 0 small 
enough and take t — 1 + s — 2 s and /3 = 1 — s. From (3.10) we obtain the following 
weighted version, where the gaining of regularity is of almost one derivative: 

Proposition 3.12. Let 1/2 < s < 1, f £ C 1_s (f2) and u be the solution of our 
problem. Then, given e > 0 it holds that u £ {LI) and 




C{SI, i 


l-s) 


Remark 3.13. The regularity estimates given in this section are sharp, in the 
sense that if we consider the problem 


f (—A) s u=linB(x 0 ,r), 
{ u = 0 in B{ 0, r) c , 


for xo £ R” and r > 0, then its solution is given by [23] 

2-2sp (n\ 

u{x) = — . „ . —-(r 2 — |a; — a;o| 2 )* in B{x$,r). 

V / r (n±2 £) r ( 1 + s ) V l 01 ) v ; 

It is straightforward to check that this function belongs to H s+ ^~ e {Ll) for all s £ (0,1), 
to Hfj^~ e e {Ll) if s £ (1/2,1) and that the parameter e can not be removed. 

3.3. Global Regularity. A direct derivation of global regularity is a simple task 
in the present context. First we present the following lemma. 

Lemma 3.14. For i < s < 1, e > 0 and u £ H s+ ^~ E {Ll), it holds 


\Du{x)\‘ 2 


n Jn<= x 


_ ^|n+2(s— i—e) 


dydx< C ( n ’ s ’ n ) \\ Du f 1 
y ~ 2s — 1 — 2e 


Proof. This is a simple consequence of the inclusion S2 C C B{x , S(x)) c for all x £ LI 
and the Hardy inequality (2.3). □ 

Combining Lemmas 3.10, 3.14, and Proposition 3.11 we have proved: 

Proposition 3.15. If 1/2 < s < 1 and f £ C^{Ll) for some (3 > 0, then the 
solution u of (2.1) belongs to H S+ ^~ E { R") for all e > 0 and 

C{Ll,s,n,/3) 

~ yfe{2s~l) 

In a similar fashion we get 

Proposition 3.16. Let 1/2 < s < 1, f £ C 1— s (f2) and u be the solution of our 
problem. Then, given £ > 0, u £ R”) and 


K/ + 2 _t £ ( 




c{n, 


l-s) 


£ 
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3.4. The case s = 1/2. Up to now, the possibility of s being equal to 1/2 
has been excluded from our analysis. In order to obtain a regularity estimate, the 
arguments to be carried are in the same spirit as before; the only issue to overcome is 
the need for /? > 0 in Theorem 3.2. In this case, the argument demands less regularity 
of the function /. Indeed, the same technique as in the proof of Lemma 3.10 gives 
u £ Ff 1_e (fl) for all e > 0, with a bound of the type 

(3-12) Mgi-e(O) < ^ £ ’ ^ II/IIl°°(Q)- 

Observe that we cannot assure u £ iJ^fl) by taking e —> 0 in the previous inequality, 
which is coherent with example (3.11). 

Moreover, in this case the space V coincides with the Lions-Magenes space H 0 q (U), 
which is strictly contained in H 1 > 2 (Q') ([31], Theorem 1.11.7). Actually, the energy 
norm is equivalent to ||m||+ ||d _1 ' 2 u||^ 2 (n). By means of the Hardy inequality 
it can be bounded by |w| ff i/ 2 +e(Q) for any e > 0. Asa consequence, FE error estimates 
for this case follow from the theory developed below for s ^ 1/2. 

4. Finite Element approximations. We restrict the analysis to FE approxi¬ 
mations of equation (2.1) to piecewise linear functions. The simpler case of Vo, which 
provides a conforming method for s < 1/2, is not addressed here in order to present 
an unified approach for the whole range 0 < s < 1. 

We assume that U = ^ where Th is an admissible triangulation of Q, made 
up of elements T of diameter hr and with px equal to the diameter of the largest ball 
contained in T. 

We require that the family of triangulations under consideration satisfies: 

(Regularity) 3a > 0 s.t. hr < apr VT £ Th, 

(Local quasi-uniformity) 3A > 0 s.t. hr < A hx' VT, T' £ Th '■ f (T T' ^ 0. 

Naturally the second condition is a consequence of the first one. In this way A can be 
expressed in terms of er. 

Consider the discrete space 

Y h ={v£Y: v\ T eVi MT&Th}- 

It is immediate to check that there exists a unique solution to the discrete problem 

(4.1) find Uh £ such that a(uh,Vh) = / fv h Vv h £ V*, 

J n 

where a is the bilinear form defined by (2.2) and that Cea’s Lemma holds in this 
context. Namely, the FE solution is the best approximation in Yh to the solution of 
problem (1.3): 

(4.2) ||u-ith||v= min ||u - t>h||v- 


It seems clear that the norm || • ||v is the ‘natural one’ for studying our problem. 
A simple trick (see Lemma 5.1) shows indeed that although it involves an integration 
over an unbounded domain, it is possible to carry the computation of the error under 
this norm by integrating in H. 
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4.1. Estimates for the Scott-Zhang interpolation operator. The next step 
towards obtaining a finite element error estimate is to find an adequate interpolation 
or projection operator. One difficult aspect dealing with fractional seminorms is that 
they are not additive with respect to the decomposition of domains. Nevertheless 
some localization is possible [19, 20]. 

Lemma 4.1. Let s £ (0,1) and f 1 a Lispchitz domain. Then, for any v £ H s { fl) 
it holds that 


lfu(n) ^ s) ^ 




y 7 ( lf dydx 
\r. - n n + 2s w 


■ rai'x'i 


where 


St := [J T'. 

T : T'nT^0 


In the following, let lihV denote the Scott-Zhang interpolator of v. Let us recall 
its basic properties [40]. 

Theorem 4.2. Let £ > 1/2, then II^ : H e (Cl) -a Yh satisfies that = Vh 

for all Vh £ Yh and 11^ preserves boundary conditions, in the sense that Hq(£ 1) is 
mapped to Y h0 := {v h €Y h : v h \ gn = 0}. 

Stability and approximability results for the Scott-Zhang interpolation in frac¬ 
tional spaces were studied in [ 12 ], where the following result is proved: 

Theorem 4.3. Given ^ < £ < 1 and 0 < t < £, then VT £ 7/, v £ H^^St) 

(4.3) \Ll h v\ H t(T) < C{n, s, a) (jiff |M| l 2 (St) + , 

and 

(4.4) Ik - fthv\\ H t(T) < C'(n,s,cr)/i^r t |u|j ? £ (ST) Vw £ H\S t ). 


The procedure towards obtaining the previous results is as follows: the stability 
of the operator (4.3) relies on basic estimates for the parametrization of T and of the 
basis functions involved; in order to obtain (4.4), it is enough to apply the stability 
estimate, the Bramble-Hilbert lemma and interpolate between some L 2 (T) and H 1 (T) 
estimates. 

The same machinery as in [12] and the previous lemma yields the following sta¬ 
bility type estimate. 

Proposition 4.4. Let T e% 

1. If s e (0,1/2], £ £ (1/2,1) and v £ H e ((l) 


|n h v(iE) - n h v(y )\ 2 


T J s T 


\x-y\ 


n+2s 


dy dx < 


<C(n,s,cr) + h? 2 s M hi(s t ) 


2. If s £ (1/2,1) and v £ H 1 ( ft) 

|n h v{x) - U h v(y)\ 2 


II 

Jt J St 


- y I 


n+2s 


dy dx < 


< C{n , s, cr) 


^T 2 S \\ v \\l 2 (St 


i 2 T 2 S W 2 


H 1 (St) 
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Remark 4.5. Regarding the behavior of the constant C in the variable s, it is 
easy to check that C ~ . 

Before obtaining error estimates for IRu we recall some facts about a well known 
key tool. Let S be an star-shaped domain w.r.t. a ball B. Introduce the polynomial 
Vk(u) of degree k with the property 

[ D a (u-V k (u)) = 0, 

Js 

for 0 < \a\ < k. In our context we need to focus on the cases k = 0,1. For instance, 
Proposition 2.2 gives at once 

Ik - V 0 {u)\\ L 2 {St) < Ch l \u\ H t(s T ), 

for 0 < £ < 1 and with a constant depending on the chunkiness parameter of St- In 
this context we can write C = C{a) (thanks to the mesh properties (Regularity) and 
(Local quasi-uniformity)). 

Since |w—7 ? o(u)|^(s T ) = M/p?(s T )j by means of the L 2 estimate and interpolation 
of spaces we obtain 

\u - Vq{u)\ H s( S t ) < Ch £ ~ s \u\ H i(g T y 

for any 0 < s < £ < 1, with a constant C = C{a). 

Similarly, using the standard Poincare inequality for functions with zero average 
together with Proposition 2.2 , we obtain for any 1 < £ < 2 

(4.5) |k - Vi{u)\\ L 2 ( St ) + hr\u - Vi(u)\ H i(s T ) < Ch l T \u\ H e^ ST ), 

with C uniformly bounded in terms of a. 

Moreover, interpolation of spaces and (4.5) give for 0 < s < 1 and 1 < £ < 2 

(4.6) \u - Vi(u)\ H s( S t ) < Ch^T s \u\ H i (s T ), 
with C bounded again in terms of a. 


4.2. Uniform Meshes. From (1) (resp. (2)) of Proposition 4.4 and approxima¬ 
tion properties of Vq (resp. Vi) for 0 < s < £ < 1 (resp. 1/2 < s < 1 and 1 < £ < 2 ) 
it follows, in an standard fashion, the local approximability inequality 

|(u - II/ l u)(a;) - {v - U h v)(y )\ 2 


II 

Jt J St 


It 

Therefore 


I 


n -\- 2 s 


dydx < C(n , s, cr)hif 


2 £— 2 s | |2 

v \H‘(St)- 


Ik - i;11v < C(n, s, a)h e s k|^ (0 )- 


Calling h = ma ~x-TeT h hr, the mesh size parameter, invoking identity (4.2), and 
combining this respectively with Proposition 3.6, estimate (3.12) and Proposition 
3.11, we have proved: 

Theorem 4.6. For the solution u of (2.1) and its FE approximation Uh given 
by (4.1) we have the a priori estimates 

||m- u h \\ v < Ve>0, if s < 1/2, 

Ik - «h||v < Ve > 0, if s = 1/2, 

Ik ~ u h \\y < ^~ E 11/II^(0) V £ > 0 , 


if s > 1 / 2 . 
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So, if h is sufficiently small, taking e = | ln/i| 1 we obtain the quasi-optimal estimates 


||m-wJv < C(s,c7)/i=|ln/i|||/|| c i_ S(n) , if s < 1/2, 

||u- u h \\y < C(a)h^\lnh\\\f\\ L ^^n), if s = 1/2, 


Hw-UfcHv < h ' ^ ln ^lll/ll^ ( o), if s > 1/2. 


4.3. Graded Meshes. The obtained approximability property of 11^ is enough 
to deal with standard fractional spaces. Nevertheless, for 1/2 < s < 1 it is possible 
to improve the convergence rate using graded meshes. It requires dealing with the 
weights already introduced in Subsection 3.2. In order to get appropriate bounds in 
these spaces we should replace the classical Poincare inequality of Proposition 2.2 by 
an improved version in fractional spaces. The term improved in this context usually 
involves weights which are powers of the distance to the boundary. On the other 
hand, in [28, Theorem 3.1] we find the following fractional improved Sobolev-Poincare 
inequality for functions with zero average 



(4.7) 


The parameters r and a can be taken arbitrarily as long as r, a € (0,1), while 1 < 
p < q < , p < n/o. In [28] the domain 12 is assumed to belong to the class of 

John domains (for a definition and properties of this class see for instance [33]); this 
class is much broader than the one of star-shaped domains. 

Now we set r = 1/2 in (4.7). For <j to be chosen, we consider p such that 
- r ff > =2 = q. Observe that this election obviously implies that p < 2, and therefore 
for all a£l, applying Holder’s inequality with exponents ^ and , 


IMUa ( n )<CI?I 2 2p , 


where 



and 


// 

Jet Jet 


q JnnB(x,^fJ) 


\x~y | 



Since for every x £ 0 and y e B{x, ^fp-) it holds that S(x, y) € ^^-,6(x) , assuming 
that er < s the second integral I 2 can be estimated as follows: 



2p(s — cr) — 2a 


dx. 
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This is finite if and only if 2p ' s 2 f p 2a > —1, and recalling the election of p we 


made, it is enough to consider 


a < 


2 n(s — a) + 2cr 
n + 2a 


Choosing a according to this restriction, we obtain that the weight in the term I\ 
must be 


2a „ 

— < 2s — 2er 
P 



Therefore, taking er = 2 (^L\) f° r some £ > 0 small enough, we obtain the Proposition 
4.7 below. The result is stated here for an star-shaped domain S. Actually, the 
constant C in (4.7) depends on the constants associated to the John domain 12. In 
the case of a star-shaped domain the John constants are easily bounded in terms of 
the chunkiness parameter. Working with a domain of diameter one, a further scaling 
argument shows the final dependence on the diameter of S. 

Proposition 4.7 (Weighted fractional Poincare inequality). Let 0 < s < 1, 
a < s and S a domain which is star-shaped w.r.t. a ball B. Then, there exists a 
constant C such that for every v € L 2 (S ) with f s v = 0, it holds 


( 4 - 8 ) IMIl 2 (S) < Cd s s 

with a constant C depending on the chunkiness parameter. 

Remark 4.8. The previous result as stated suffices to fulfill our needs. Never¬ 
theless, from the standard version with s = 1 one might expect (4.8) to hold even if 
a = s. This is indeed the case, however we were unable to produce a proof as short 
as the one given here for a < s. 

Now we want to exploit the previous proposition together with the a priori esti¬ 
mate of Proposition 3.12. Since the weights under consideration vanish only on the 
boundary of the domain we need to rely on (4.8) just for patches St touching 912. 
Actually, for them we obtain the following improved version of (4.5), derived using 
Proposition 4.7 instead of Proposition 2.2 

II U - 'Pi('u)|| i 2( 5T ) + h T \u - Vi(u)\ h -l^St) ^ a \ u \Hi(S T )’ 

where 1 < i < 2 and a < i — 1. Taking 1/2 < s < 1, £ = 1 + s — 2e, and a = 1/2 — e 
we obtain the analogous of (4.6) 


■ ~ T-’i{'u)\h*(s t ) ^ CTij / 2 £ \u\ h i+b~ 


1/2 -"(St)- 


In particular, this property of Vi and the stability estimates (see (b) in Proposition 
4.4) yield 


|(u - n /l i/)(a;) - (v - U h v)(y)\ 2 


T J S T 


- y I 


n+2s 


dydx<C{n,s,a)h T e \v\ ± +a - 2 e (q r 

n i/2-e ) 


This approximability property is particularly useful for patches St touching the 
boundary of 12. For these it must be recalled that dist{x,dST) < dist(x,d 12) for 
x £ St- 

The following is standard (see [25, Section 8.4]). We assume that, in addition 
to (Regularity) and (Local quasi-uniformity) our meshes enjoy some extra properties, 
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denoted below with ( H ). First, we pick an arbitrary mesh size parameter 0 < h and 
define, for e small enough, a number 1 < /x = 2/(1 + 2er) < 2. 

Property (H): assume that for any T £ % 

• If T fl 512 ^ 0, then hr < Ch M 

• Otherwise hr < C/i disi(T, 512) b* -1 )/^ 

Using the estimates for 11/! given in Subsection 4.2 when St 0 512 = 0 and the a 
priori estimate Proposition 3.12, we can conclude, for graded meshes obeying (H) 
and 1/2 < s < 1, that 

II u ~ U /Jv < 2s—1 ^ ^\/l l n ^-lll/llci- s (^)- 

If the mesh parameter h can be appropriately related to the number N of nodes of 
the mesh then it is possible to obtain quasi-optimal order of convergence. 

Theorem 4.9. Let 1/2 < s < 1 and assume that the FE triangulation 7h satisfies 
conditions (Regularity), (Local quasi-uniformity) as well as the grading hypotheses 
( H ). If the mesh parameter h behaves like h ~ N \/ n , N being the number of mesh 
nodes, then for the solution u of (2.1) and its FE approximation Uh given by (4.1) 
we have the following a quasi-optimal a priori estimate 

ll« - «*llv < ^J^Y^iV- 1 /n yfhOVf||/|| cl - 8 ( Q ). 


Remark 4.10. In the next section we show a concrete 2D example in which 
meshes of the kind required in previous theorem are constructed. 

5. Implementation details and results. Numerical computation of solutions 
of (1.3) has as main difficulties the fact that a singular kernel is involved, and that 
integrals over the whole R ra must be calculated. 

Now we will comment some features of the implementation, more details can be 
found in [1], Let be the nodal basis of Y h- Due to the linearity of the fractional 
Laplacian, we just need to solve a system KU = F, where the right hand side vector 
F = (/,;) is assembled straightforwardly because 


fi 



dx. 


The challenging task is to accurately compute the stiffness matrix K = ( Kij ), given 

by 


Kij = 


(<Pi(x) - l Pi(y))( l Pj{x) - Vj{y)) 


- y I 


n+2s 


dx dy. 


Splitting Q = (12 x 12)U(12 x I2 c )U(r2 c x 12) and taking into account that the interactions 
in 12 x f 2 c and 12 c x 12 are symmetric respect to x and y, we get 


K tj = 


Q,xQ 


+ 2 


(¥>i(x) - <Pi(y))(<Pj(x) - ipj(y)) 
\x - y\ n + 2s 

ipi(x)(pj(x) 


dx dy - 


f 2 xO c 


\x - y I 


n+2s 


dx dy. 


By making a double loop over the elements of the triangulation, the integrals 
above can be computed. The quadrature rules employed for computing the integrals 
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over two elements T and T' (with the possibility that T = T') are analogous to the 
ones presented in Chapter 5 of [39]. The advantage of applying the transformations 
presented in that book for this problem is that they convert an integral over the prod¬ 
uct of two elements into an integral over [0,1] 4 , in which variables can be separated 
and the singular part can be solved analytically. The integral involving Cl x Cl c is 
computed resorting to integration in polar coordinates, taking into account that it 
will be nonzero only if supp(</?i) (T supp(^Jj) ^ 0. 

5.1. Numerical Results for Uniform Meshes. In first place, numerical so¬ 
lutions of problem (3.11) were obtained for n = 2, xo = 0 and r = 1, and for several 
values of s. The computation of the error in the V norm is easily achieved by using 
the following. 

Lemma 5.1. It holds 

||u - u/JIv = f{x)(u(x) - u h {x)) dx 

Proof. It is an immediate consequence of the orthogonality condition 
a(v h ,u-u h )=0 Vr/, £ Vft. 

Indeed, from it we obtain 

||«- u h \\y = a(u- u h ,u- u h ) = a(u, u — u h ), 

and the equality follows by (2.1) and (2.2). □ 

Although the computation involved in this lemma is subtle in general, in this 
particular case it can be carried out exactly since / = 1 on and a closed formula 
for J n u is easy to get while the exact value of f n Uh can be numerically evaluated. 

Several orders are shown in Table 5.1; these results are in accordance with the 
estimates in Theorem 4.6. In Figure 5.1 computational errors for s = 0.5 and s = 0.7 
are shown. 



Value of s 

Order (in h ) 

0.1 

0.497 

0.2 

0.496 

0.3 

0.498 

0.4 

0.500 

0.5 

0.501 

0.6 

0.505 

0.7 

0.504 

0.8 

0.503 

0.9 

0.532 


Table 1 

(Uniform Meshes) Computational rates of convergence for problem (3.11), measured in the 
norm || ■ ||y. The mesh parameter is the actual size of the elements. 


As a second example, take s > 1/2 and consider problem (1.3) posed on the 
interval Cl = (—1,1), with exact solution u(x) = sin(7rx)x(_ 1 ,!)(*), namely: 


(5.1) 


f (—A) s u = (—A) s sin(7ra:) in (—1,1) 

( u = 0 in (—oo, — 1) U (1, oo). 


Since the solution for this problem is smooth in (—1,1), the convergence in the energy 
norm would be expected to be of order 2 — s. Some results are shown in Table 2, 
where it can be seen that these orders are indeed achieved. 
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Fig. 1. Computational results for problem (3.11) using uniform meshes. The left panel shows 
the rate for s = 0.5 and the right one for s = 0.7. In both cases, the rate is ~ 0.5, as predicted by 
Theorem f.6. 


Value of s 

Order (in h) 

0.6 

1.4028 

0.7 

1.2993 

0.8 

1.2002 

0.9 

1.1002 


Table 2 

Rates of convergence for uniform meshes in the norm || • ||y for problem (5.1) and s > 1/2. 


5.2. Numerical Results for Graded Meshes. For the same 2D problem of 
the first example we show how to build appropriate graded meshes. Our domain Q is 
the unitary disk. Therefore, we may pick a positive integer M and define an increasing 
sequence of radii r» := 1 — (l — -jjV for 1 < * < M. We can mesh the complete disk 
0 by meshing each subdomain flj = {x £ 0 : ?y_i < \x\ < Ty} with uniform elements 
of size h,T = hi = Ti — r;_i (see Figure 5.2). proceeding in that way it is possible to 
compute the final number of nodes N ~ y~) ._ 1 1/hi. It is a simple exercise to check 
that if fi < 2 then N ~ M 2 . The previous construction ensures that conditions 
(Regularity), (Local quasi-uniformity) and hypotheses ( H ) hold, taking h = 1/M. 

Table 3 shows numerical results for this case. The accuracy is in full agreement 
with that predicted in Theorem 4.9. 



Fig. 2. Left: graded mesh with M = 15 and fi = 2 — e. Right: uniform mesh with M = 15 and 
M = 1- 


Remark 5.2. Taking into account the restrictions (3.9), it is possible to achieve 
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Value of s 

Order (in h) 

0.5 

1.066 

0.6 

1.051 

0.7 

0.990 

0.8 

0.985 

0.9 

0.977 


Table 3 

(Graded Meshes) Rates of convergence in the norm || ■ ||y for problem (3.11) and s > 1/2. The 
mesh parameter h behaves like N ~ 1 / 2 , N being the number of nodes. 


differentiability orders between 1/2 + s < £ < 2 by choosing adequate weights. At this 
point, the reader might ask whether the order of convergence (with respect to N) could 
be improved by considering a different value of £ and following the grading approach 
presented at the beginning of this subsection. This is not the case; actually the choice 
we made yields the best possible order w.r.t. the number of nodes with minimum 
grading requirements on the mesh. 

Indeed, it is simple to check that, for a given regularity I, the optimal choice for 
the grading parameter is p = 2(1— s). Recall identity (3.10) and the results shown in 
Section f, which give 

< Ch e ~ s \u\ H e a{n) , 

where h is the mesh parameter. 

If we restrict to I < 1 + s, then p < 2 and the number of nodes is N ~ /i -2 . 
Therefore, as £ increases there is a gain of order without an increment in the total 
number of nodes and the error behaves like . Within this range, the choice 

1=1 + s — e is optimal. 

On the other hand, if we consider £ > 1 + s then p > 2 and it is simple to check 
that in this case N ~ h~^. Here the gain of order one might expect due to the increase 
in differentiability is compensated by the cost of having to increase the weight power, 
which implies a growth in the number of nodes. In the whole range £ £ (1 + s,2) we 
obtain that the error behaves like N~ 1 ^ 2 In N. 

6. Conclusion. In this paper, a complete Finite Element study of a fractional 
Laplace equation is carried out. First it is shown that recent Holder regularity results 
for this problem [38] can be used to provide a priori estimates in weighted fractional 
Sobolev spaces, within which the FE setting can be straightforwardly adapted. In 
particular, some of these estimates measure in a precise way the singular behavior of 
solutions near the boundary. Borrowing techniques from the BEM it was found that 
the singular kernel arising in this problem can be accurately handled. The FE method 
is implemented in one and two dimensions, where uniform as well as tailored graded 
meshes are proposed. Error estimates for the Scott-Zhang interpolation operator 
in fractional weighted spaces are obtained by introducing an appropriate version of 
the improved Poincare inequality. These error estimates are used to prove optimal 
order of convergence of the FE method in the weighted fractional context. Numerical 
experiments are presented delivering orders of convergence in full agreement with our 
theoretical predictions. 
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